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Abstract 

The Generalized Method of Cells (GMC),a micromechanics based 
constitutive model, is implemented into the finite element code MARC 
using the user subroutine HYPELA. Comparisons in terms of trans- 
verse deformation response, micro stress and strain distributions, and 
required CPU time are presented for GMC and finite element models 
of a fiber/matrix unit cell. GMC is shown to provide comparable pre- 
dictions of the composite behavior and requires signficantly less CPU 
time as compared to a finite element analysis of the unit cell. Details 
as to the organization of the HYPELA code are provided with the 
actual HYPELA code included in the appendix. 

1.0 Introduction 

In the design of aerospace components, the use of composite materials 
such as metal matrix composites (MMC) are being investigated as a means 
of obtaining increased operating loads, improved durability, and decreased 
weight. For example, MMC’s are being considered for critical engine compo- 
nents, such as rotors, blades and nozzles which are subjected to high temper- 
ature environments under complex load histories. From the analysis stand- 
point, the use of nonlinear finite element methods will be required for the 


1 


design and structural response predictions of these components. In addition, 
there is increased interest in developing life prediction methods to predict 
the service life of these components which will also require accurate finite 
element analysis. 

When analyzing a structure using the finite element method a key ingre- 
dient is the constitutive model as it determines how accurate the global re- 
sponse of the structure will be predicted. For structures composed of isotropic 
materials, there is considerable experience with plastic and viscoplastic con- 
stitutive models and their implementation into the nonlinear finite element 
solution. On the other hand, the development of effective constitutive mod- 
els for the nonlinear material response of a composite is relatively new. The 
primary goal in formulating a constitutive model for composite materials is 
the development of the relationship between the macro strain and the macro 
stress, i.e. the effective macro stiffness. There are two approaches [2] in mod- 
eling the composite material: 1) a macro, equivalent, homogeneous material 
model, 2) a micromechanics based model. Considering the operating envi- 
ronment that the components will be subjected to, using a micromechamcs 
based material model may be necessary in order to include the many non- 
linear mechanisms that occur within the composite, such as, fiber/matrix 
debonding, localized matrix yielding, matrix cracking, oxidation effects, etc.. 
These type of effects are believed to be especially important for the life 
analysis of the component. There has been debate whether or not using a 
micromechanics based constitutive model in the finite element method is fea- 
sible with regards to computational expense. This point will be addressed in 
the following section by examining some of the various approaches currently 
available. 

2.0 Computational Aspects 

The general solution algorithm for a nonlinear finite element analysis 
will be discussed in terms of the computational effort required at both the 
global (structural) level and the local (element) level. This discussion will be 
used to show how the use of micromechanics based models may or may not 
significantly increase the computational expense of a nonlinear finite element 
analysis. 
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(1) 


The nonlinear finite element equilibrium equation takes the form, 

[K] r {AU} = {AR} 

where [K] r is the structural tangent stiffness matrix and the term {AR} is 
the out-of-balance force vector, i.e., 

{AR} = {F} - p l'lB\T{<r},\3\^dt,drds (2) 

Global equilibrium is satisfied by iterating on eqn. (1) until the out-of- 
balance force vector { AR} is effectively zero (within some predefined conver- 
gence tolerance). This global equilibrium equation therefore constitutes the 
first level of iterations. 

From the out-of-balance vector {AR}, eqn. (2), the integral is evaluated 
at each integration point in the structure. Specifically, at each integration 
point, the global finite element solution provides the current strain state and 
the current element stress u is obtained from the evaluation of the local 
constitutive model. That is, the constitutive model calculates the current 
stress and should provide the current material stiffness at that point. Thus a 
micromechanics model needs to be able to take the macro/composite strain 
and, after performing a local stress/strain analysis on the RVE, provide the 
current macro/composite stress state and equivalent composite material stiff- 
ness. Having the current equivalent composite material stiffness will allow 
the use of a tangent stiffness matrix in the global equilibrium equation (1). 
The advantage of using a tangent stiffness is the increased rate of conver- 
gence in the global equilibrium iterations, with a quadratic rate obtainable 
if a properly formulated local tangent stiffness is utilized. 

It is apparent that the next important level of computations is at the 
integration point. Here, one may or may not encounter local level iterations 
which will depend on the type of constitutive model being used. For clarity 
in the following discussion, assume that the material at the integration point 
in question is isotropic. 

If for example, a classical plasticity model is used, the radial return al- 
gorithm, which is commonly used, requires local level iterations in order to 
return the stress to the yield surface, an example of which is given in [11]. 
Another class of material models would be the continuum-based viscoplastic 
models. These rate dependent models contain flow and evolution equations 
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which must be integrated in time. This integration is performed using ei- 
ther an explicit or implicit integration method. If an explicit integrator is 
used, for example forward Euler, no additional local iterations are required. 
But it is well documented that such explicit methods are conditionally sta- 
ble and may require small global time steps in order to achieve convergence. 
Solving a highly nonlinear structural problem with many small time steps 
becomes computationally expensive. As a result, the use of implicit integra- 
tors are increasingly being used since they are far more robust for integrating 
the numerically “stiff” viscoplastic rate equations and provide unconditional 
stability, thus allowing larger time steps to be used. But, due to their implicit 
formulation, local interations are required. 

From the above discussion, it is apparent that two levels of iterations are 
typically required to solve highly nonlinear finite element problems: 1) global 
equilibrium, and 2) local constitutive model iterations. These multiple levels 
of iterations are what makes nonlinear finite element solutions inherently 
computationally expensive. Therefore, it is critical that any micromechanics 
constitutive model be accurate yet efficient. 

Any micromechanics method which is based upon a finite element RVE 
becomes in effect a finite element analysis within a finite element analysis. 
A complete nonlinear finite element analysis must be performed to take the 
current macro strain state provided by the finite element solution and pro- 
duce the current stress state in the RVE. Note that at this point, only the 
micro stress/strain states exist. Additional computations are still required 
to convert the micro (constituent) stresses to an equivalent macro compos- 
ite stress state and calculate the equivalent macro material stiffness. An 
example of a finite element based method is that presented by Wu et. al. 
[15]. In this method a RVE for a periodic hexagonal array (PH A) finite el- 
ement model is used to preform the micromechanics analysis and calculate 
the instantaneous macro composite material stiffness. The PHA model in 
effect serves as the local constitutive model and is called at each integration 
point in the structure. The ultimate objective of a micromechanics analysis 
is the calculation of the equivalent macro material stiffness. In this regard, 
a so-called equivalent homogeneous volume, EHV, is introduced. The EHV 
is defined by the equivalent macro instantaneous material parameters which 
are considered as unknowns. These material parameters are determined by 
applying the same periodic boundary conditions to the EHV as are applied 
to the RVE. Iterations are then performed until equilibrium is attained such 
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that the EHV and RVE have the same total energy change. One can see 
that the level of computational effort is increased. In addition, each finite 
element in this RVE has a constitutive model which must be evaluated and, 
as discussed above, may or may not require local iterations in its calculations. 
It is stated in [15] that accurate calculation at each material point requires 
substantial effort which would need to be repeated many times in a typical 
nonlinear finite element analysis. 

Another, more general, finite element based approach is that of Fish and 
Wagiman [7] which is applicable to both periodic and non-periodic domains. 
That is, the solution or the microstructure may be periodic or non-periodic. 
In this approach, the homogenized material properties are based upon the 
displacement field of the unit cell which is obtained by solving the finite 
element model of the unit cell for the six different loading conditions. In ad- 
dition, in this “multiscale” finite element method, two finite element meshes 
may be used. A “macro” mesh is used for the global problem which is the 
part of the structure where the composite microstructure and the solution is 
periodic. A second “micro” mesh is used for the local problem where either 
the solution or the microstructure is not periodic. This produces a coupled 
global/local problem which, as described in [7], requires an iterative solution 
process. 

From the above discussion, it is apparent that while the use of finite 
element based micromechanics models provide significant capabilities, the 
associated computational expense (time and memory) and complicated algo- 
rithms make them less attractive. The use of parallel computers and paral- 
lelized/optimized code [13] in order to decrease computation time limits the 
usefulness of the method since such sophisticated computer resources are not 
commonly available. 

In order to reduce the computational expense, an efficient micromechan- 
ics model should be more analytical in nature in which closed form expres- 
sions are available for the macro/micro relationships. This class of models 
are based upon micromechanics but include various assumptions in order to 
derive semi-analytical expressions for the effective composite stiffness and 
relationships between macro and micro stress and strain. 

Svobodnik et. al. [12] use the vanishing fiber diameter model of Dvorak 
and Bahei-El-Din [6] in the context of the finite element method. In the 
vanishing fiber model, the fibers are assumed to behave elastically to fail- 
ure and any effects of fiber interaction are neglected. The result is a simple 
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model which gives the equivalent macro response of the composite expressed 
in terms of the constituent properties and their volume fractions. Svobodmk 
[12] demonstrates that the model may be implemented into the ABAQUS 
finite element code with relative ease. There are limitations of the vanishing 
fiber model. The first is neglecting any fiber interaction effects. It has been 
shown that packing arrangement, i.e. square, hexagonal, etc., does have an 
effect on the transverse response of the composite [3]. A second limitation of 
the model is the inability to capture specific micro stress/ strain distribution 
effects such as localized yielding in the matrix [12] due to its assumption of 
homogeneous stress and strain fields in the fiber and matrix. Such a lim- 
itation is important if accurate micro stress and strain is necessary as in 
the case of damage /life analysis of the composite. With these limitations, 
it is believed that any benefit of using a micromechanics based model have 
been lost. The choice of micromechanics-based models is either a simplified 
model which provides the necessary macro composite stress strain relation 
but because of its simplified form is not able to provide micro stress strain 
information in sufficient detail. On the other hand, finite element based mi- 
cromechanics models do provide detailed (depending on mesh density) fiber 
and matrix stress and strain distributions but at the cost of significant com- 
putational expense in terms of the microlevel analysis and the calculation of 
the current macro composite material stiffness. Thus the “best” microme- 
chanics model would be a combination of the above two approaches, that is, 
it would provide sufficiently detailed microlevel stress and strain distributions 
and a macro composite material stiffness in an efficient manner. 

In this regard, the Generalized Method of Cells (GMC) micromechanics 
model of Paley and Aboudi [9] is believed to satisfy the above requirements. 
GMC is an extension of Aboudi’s earlier method of cells [1] in that GMC 
now allows for more than 4 subcells to describe the fiber and matrix phases. 
By having a variable number of subcells provides significant flexibility in 
modeling a variety of RVE’s, for example different fiber shapes, fiber packing 
arrangements, inclusion of a fiber /matrix interface layer, etc. Some typical 
RVE’s constructed using GMC are shown in Figure 1. In addition, one of 
the noteworthy features of the GMC model is that it provides an analytical 
constitutive expression for the composite as well as macro-micro relations 
for stress and strain. Thus no iterations are required, as in a finite element 
based model, to obtain the current macro composite material stiffness. Yet, 
as will be shown in the following sections, GMC provides detailed microlevel 
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stress and strain distributions which are comparable to those obtained from 


a finite element RVE analysis. 

Some details of GMC pertaining to its implementation in the HYPELA 
subroutine in MARC will be briefly outlined below. Additional theoretical 
details may be found in [9]. 

3.0 Generalized Method of Cells: The HYPELA Subroutine 

As mentioned in section 1 , one of the uses for an efficient micromechanics 
model would be in the finite element method. In order to access the re- 
quirements of a micromechanics model in terms of CPU time and overhead, 
i.e. memory/storage requirements, the GMC model has been implemented 
into the finite element code MARC through the use of the user subroutine 
HYPELA [8]. The subroutine HYPELA is provided for users to implement 
their own constitutive models into MARC. The implementation of GMC into 
a HYPELA format was relatively straightforward. The main subroutine is 
HYPELA which contains all of the micro and macro relationships that com- 
prise GMC. A branch is made to the desired constitutive model for each 
subcell material. The equations below are presented in the same order in 
which they appear in the HYPELA subroutine. One of the key parts of 
GMC is the formation of the A-matrix. This matrix consists of submatrices 
which contain quantities defining the material constants and the geometry 
of the individual subcells. The A-matrix is partitioned in the HYPELA 
subroutine as follows, 


A M 0 Am 
A q J 0 


(3) 


Let the submatrix A* be defined as, 



(4) 


where A* is a square matrix with dimensions NpX N y . Also define K as, 



( 5 ) 
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and A p as 


( 6 ) 



In the above matrices, the matrix A m involves the elastic material properties 
of the subcell and the matrix A G involves the geometry of the subcell and the 
matrix J relates the average composite strain to the average subcell strain. 

The only necessary modification to the code included in Appendix A is 
the linking of a matrix routine for the inversion of the A matrix. It will be 
noted that in the present code, two Unpack subroutines, i.e. SGEFA and 
SGESL, were originally used for this purpose. The calls to these subroutines 
are left in the code to show the proper location for the matrix inversion 
subroutine call. Also note that since only the elastic material constants are 
being used, it is only necessary to form the A matrix at the initial time step. 
On the other hand, if viscoplasticity (tangent stiffness) or damage effects are 
to be included in the stiffness, the A matrix would need to be calculated at 
each time step. 

Once the A-matrix has been formed, the effective macro composite elastic 
stiffness may be calculated from, 

C = Tj E E A«*> P) 

^ p=l 7=1 

in which is the subcell elastic stiffness. Note that calculating the crit- 

ical macro stress/strain relationship necessary for the global finite element 
solution only required straightforward matrix multipHcation operations. 

Another advantage of GMC is that the inelastic response of the con- 
stituents may be modeled using either a plastic or viscoplastic constitutive 
model. The quantity is the inelastic strain rate in the subcell and cal- 

culated by calling the appropriate constituent constitutive model subroutine. 
In the current implementation, an isothermal Bodner-Partom [5] viscoplas- 
tic material model is included. Additional material models can be included 
with relative ease. Note that both the viscoplastic material constants for 
the Bodner-Partom model and the elastic material constants are defined in 
the HYPELA/GMC subroutine. The viscoplastic constants are stored in the 
two dimensional array VPR.OB and currently is dimensioned for two differ- 
ent materials having a maximum of 12 viscoplastic constants per material. 
Note that the constitutive model calculates the inelastic strain rate for each 
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subcell. That is the matrix subcell inelastic strain, is calculated from 

the flow law, 





( 8 ) 


From the individual subcell inelastic strain rates, the macro composite in- 
elastic strain rate, e , is calculated from, 


a 

e 


= C- 1 - E E A p( ^>cf - 

hi 


0=1 7=1 


( 9 ) 


and is stored in the array EPG. 

The average subcell total strain rate, is stored in the array ES and 

is calculated from the relation, 

^7) = A -i(07) K | + A p(/37) el (10) 


The average subcell stress rate, is stored in the array SS and is found 

from = OiiJ 7 ’ - (11) 

where C\^i is the elastic material stiffness of the (Pj) subcell. 

The macro composite strain rate, e is stored in the array EG and is 
calculated from the solution of the global (structural) finite element solution. 
Note that in HYPELA, MARC provides the current strain increment in the 
array DE and the rate is simply obtained by dividing each component by At. 
Finally, the macro composite stress, b, is obtained from the expression, 

b = C(E-e) ( 12 ) 


and is stored in the array SG. 

Before leaving the subroutine HYPELA it is necessary to integrate the 
above rate equations. Presently, the calculations in HYPELA are performed 
in a sequence that is consistent with a simple forward Euler integration 
scheme. The current solution at time t + At is found from, 

X t+At _ x t _|_ Ati| t (13) 
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where x is an array containing all of the macro and micro variables and the 
array x contains the rates which are evaluated based upon the previous con- 
verged solution at time t. In the subroutine HYPELA the array SIGAL is 
used to store the current total values of the stress , strain and internal state 
field variables, while the array DSIGAL stores the rate of the stress, strain 
and internal state field variables. These arrays are arranged identically. The 
first 37 positions are used to store the macro (composite) state variables 
which consist of total strain (6), Cauchy stress (6), inelastic strain (6), room 
for two internal state variables (12), thermal strains (6) and current tem- 
perature, in this order. The micro (sub cell) quantities are then stored in 
’’blocks” of 36 which are organized according to the description given above 
for the macro quantities with the exception that no temperature variable is 
stored. Presently, the arrays are dimensioned for a maximum of 49 subcells. 

When the user subroutine HYPELA is used, MARC solves for global 
equilibrium using the initial strain form of the equilibrium equation, i.e., 

[K e ] {AU} (i) = 

{F} - p j\ J\ {a}i\3\-^dtidrds - ^[B] 3 'IC]*{Ae} / <iF (14) 

where [K e ] is the global elastic stiffness. As a result, the macro composite 
stress is placed in the vector S, the current inelastic strain increment is placed 
in the vector G, and the macro composite elastic stiffness is placed in the 
array D so that they may be passed back to the main MARC program. 

4.0 Generalized Method of Cells Evaluation 

In order to access the accuracy of GMC, a series of comparisons were 
performed with finite element meshes representing the fiber and matrix unit 
cell. The finite element meshes used are shown in Figure 2. Material con- 
stants for a SiC fiber and Ti-15-3 matrix at 426C were used in the present 
analysis. The unit cell was subjected to a tensile load under strain control 
with the strain rate being 0.001/s to a maximum macro transverse strain of 
0 . 02 . 

4.1 Macrolevel: Deformation Response 
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The results for longitudinal response were identical for all of the finite 
element meshes and the GMC models. This is to be expected since even 
a simple rule of mixtures model can adequately predict the longitudinal re- 
sponse. However, the transverse response requires more accurate modeling 
of the unit cell. The transverse deformation response for the composite is 
shown in Figure 3. 

From Fig. 3, the 4 finite element mesh grossly over-predicts by approxi- 
mately 42% the transverse response as compared to the 4-cell GMC model. 
The response for the 16 finite element mesh is approximately 4% too stiff, 
while the response for the 64 finite element mesh is approximately 2% softer. 
Note however that the 49-cell GMC model gives a softer response than the 
4-cell model by approximately 8%. As a result, a 196-cell GMC model was 
analyzed and its response matched that of the 49-cell model. Thus the 49-cell 
model has converged to the correct transverse response of the composite. The 
272 finite element mesh is approximately 5% stiffer than the 49-cell GMC 
model. Finally, a finite element mesh of 1088 elements was used to match 
the response of the 49-cell GMC model. A note should be made that it may 
be possible to match the 49-cell response using a smaller finite element mesh. 
As may be noted from figure 2, the 1088 mesh was obtained by simply re- 
fining the 272 mesh by subdividing each element into four elements. Thus 
the 1088 element mesh may not be the “optimal” mesh using the minimum 
number of elements required to converge with the 49-cell GMC model. In 
any case, a finite element mesh larger than 272 elements is required to match 
the response of a 49-cell GMC model. The implications of requiring such a 
large finite element model will be discussed in section 4.3. 

4.2 Microlevel: Stress and Strain 

The other aspect of choosing a micromechanics model is how accurate 
the microlevel fields are predicted. Of particular interest are the stress and 
inelastic strain distributions. Figure 4 shows the effective stress, distri- 
bution for the 49-cell model and the 64 finite element mesh. Note how the 
distributions compare in both distribution and magnitude. Figure 5 shows 
the effective inelastic strain distribution for the GMC and finite element mod- 
els. Again, the distribution and magnitudes are comparable. One interesting 
note is with regards to the in-plane shear stress distribution. Inherent in the 
GMC method is the uncoupling of the shear and normal stress components. 
That is, for a given applied normal stress, a shear stress will not be produced. 
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As a result, there is a state of zero shear stress in the unit cell for the present 
tensile test. On the other hand, there is a shear stress present in the finite 
element solution. Figure 6 shows the shear stress distribution for the finite 
element model. It is interesting to note that even with the absence of any 
shear stress component, GMC still accurately predicts the macro transverse 
deformation response. This maybe attributed to the fact that the effective 
stress distribution does match that of the finite element solution [4]. Since 
the viscoplastic model is “driven” by J 2 , the deformation response should be 
comparable. However, the implications of GMC’s lack of shear stress distri- 
bution may surface when considering other aspects such as damage modeling. 
For example, fiber /matrix debonding which may be a function of the local 
shear stress in the interface region. Such issues are currently being addressed 
in ongoing research. 

4.3 Computation Time Comparisons 

The goal is to be able to analyze composite structures at both the macro 
and micro levels [2]. That is, to analyze the response of the structure 
(macrolevel) and then, at various points in the structure, go to the microlevel 
and determine the stress, strain and other internal state variables in the con- 
stituents of the composite. This type of macro/micro analysis capability is 
believed to be important when trying to capture the damage mechanisms oc- 
curring at the fiber/matrix level which in turn affect the life of the structure. 

In the context of the finite element method, these micro analyses would 
be performed at each integration point. Recall that when forming the fi- 
nite element local stiffness matrix, a constitutive model is called at each 
integration point where for a given strain state, the corresponding updated 
stress state and the material stiffness is returned. Thus, the 49-cell GMC 
model, contained in HYPELA, would be evaluated at each integration point. 
Similarly, one of the finite element models, either the 272 or 1088 element 
mesh depending on the desired accuracy, would have to be evaluated at each 
integration point. 

To study GMC, only a single finite element was used. Therefore, the total 
CPU time was then divided by eight (eight integration points per element) 
to obtain the approximate CPU time at one integration point. On the other 
hand, since the finite element model of the fiber/matrix unit cell would have 
to be analyzed at each integration point, the total CPU time was used for 
comparison to GMC. 
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With the above in mind, direct comparisons in terms of CPU time with 
the 4, 16, 64, 272 and 1088 finite element models of the fiber /matrix unit 
cell and the 4 and 49 cell GMC models were performed. A second reason 
for using a 49-cell model was that if other RVE’s are used to represent other 
fiber architectures (hexagonal pack, interface layer, etc.) the GMC RVE 
would consist of 49 cells on average [3]. 

All of the results presented in the Tables were generated on the Cray- 
YMP at NASA Lewis Research Center. Also, all of the results presented have 
been normalized with respect to the CPU time for the GMC model under 
comparison. In Table 1, if the 49-cell GMC model requires 1 CPU unit, 
a finite element analysis using 1088 elements that provides the equivalent 
accuracy, in terms of transverse displacement, requires approximately 3550 
CPU units. This clearly demonstrates the penalty one has to pay in terms 
of the overhead required for the finite element method and the substantial 
savings in CPU time the GMC method provides. One must also realize that 
the solution times reported here are only for the stress analysis portion of 
the unit cell. Determining the current equivalent inelastic macro and micro 
properties has not been performed. This phase will require what is believed 
to be significant computation time. Further details of the calculations that 
would be required in this phase may be found in [14]. 

An alternative approach to the above micromechanics models are the 
phenomenological continuum-based models. This class of models treat the 
composite as an homogenous anisotropic material in which no distinction is 
made between the fiber and matrix constituents. Therefore no micro stress- 
strain information is available. One of the arguments for using such models 
is that they provide the necessary macro relationships required in a struc- 
tural analysis and are believed to require the minimum amount of compu- 
tation time. For comparison purposes, Robinson’s viscoplastic, transversely 
isotropic continuum-based model [10] was characterized for the SiC/Ti-15- 
3 composite. The viscoplastic model was then implemented into MARC 
through the HYPELA subroutine. Again as in the case of the GMC model, 
a single finite element was used and the total CPU time was divided by eight 
to obtain the approximate CPU time for one integration point. The same 
tensile strain control test, as described above, was performed. Comparisons 
of the CPU time for the GMC, finite element and continuum-based models 
are given in Table 2. Note that the continuum model is 20% faster than the 
4-cell GMC model and 80% faster than the 49-cell GMC model. Of course, 
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one needs to keep in mind that the GMC models provide additional infor- 
mation on the micro (constituent) level that the continuum-based models do 
not provide. 

5.0 Conclusions 

The Generalized Method of Cells micromechanics model has been im- 
plemented into the finite element code MARC through the use of the user 
subroutine HYPELA. The implementation was found to be rather straight- 
forward requiring no more effort than other implementations of “standard 
constitutive models, for example the continuum-based viscoplastic model of 
Robinson which was also used in a HYPELA subroutine for this work. 

It was found that the finite element method requires a significant num- 
ber of elements in order to adequately predict the transverse deformation 
behavior of the composite. As a result, the computation time far exceeds 
that required by the GMC model. An important point is that the CPU 
time reported for the GMC model includes the calculations for the neces- 
sary macro stress-strain relations, i.e. current RVE material stiffness. While 
CPU times reported for the present finite element analyses are only for per- 
forming the nonlinear stress analysis on the RVE. Additional calculations, 
producing larger CPU times, are required to obtain the current RVE mate- 
rial stiffness. The speed of GMC is attributed to the analytical expressions 
for the macro /micro stress strain relationships and for the macro composite 
stiffness. 

The GMC model has been shown to provide accurate effective stress, ( J-i 
form), stress and inelastic strain distributions. One important observation 
is the lack of any micro shear stress components for the GMC model. As 
noted, this does not appear to affect the deformation response predictions of 
GMC, but how this limitation would affect particular life prediction damage 
models, which may depend on specific micro stress components, needs to be 
investigated further. For example, how would an interfacial debonding crite- 
ria, which may be a function of the local in-plane shear stress, be affected? 
Research in the area of life analysis is currently in progress. 
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Table 1: CPU Time Comparisons 


GMC 

Finite Element 

case 

CPU time 

mesh size — > 

4 

16 

64 

272 

1088 

4-cell 

1.(1. 17s) 

CPU time 

21 

25 

115 

1130 

15000 

49-cell 

1. (5.07 s) 

CPU time 

5 

6 

27 

261 i 

3550 


Table 2: Model Comparisons 


case 

GMC 

F.E. (272) 

Continuum 

4-cell 

1. (1.17s) 

1130 

0.8 

49-cell 

1. (5.07s) 

261 

0.2 
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Figure 1: GMC Representative Volume Elements 
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Figure 2: Unit Cell Finite Element Meshes 
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Figure 3: Transverse Deformation Comparison 
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Figure 4: Effective Stress Distribution 
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Figure 5: Effective Inelastic Strain Distribution 
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Appendix A: 

HYPELA subroutine for GMC 


C ############################################################################ 
SUBROUTINE HYPELA (D , G , E , DE , S , TEMP , DTEMP , NDUMM , N , NN , KC , MAT , NDM , 

& NDUMMM) 


C 

IMPLICIT REAL*8 (A-H.O-Z) 

C 

include ’ /tpsw/marc/ common/ elmcom ’ 
include ’ /tpsw/marc/common/concom ’ 

C 

c \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\A/\/\AAAAA 
c 

C NBM - NUMBER OF SUBCELLS IN 2-DIRECTION 

C NBG - NUMBER OF SUBCELLS IN 3-DIRECTION 

C NMM - TOTAL NUMBER OF DIFFERENT MATERIALS 
C 

PARAMETER (NBM=7 , NGM=7 , NMM=2) 

PARAMETER(NSM=NBM*NGM,NAM=6*NSM,NAM12=2*6*NSM+6,NMM36=36*NMM) 
PARAMETER(SQ2=1. 414215, PAI=3. 141592) 


C 

LOGICAL FORDIS , HEAT , THME .NINDPT , PATRAN 


C 


COMMON /MLGIC/ FORDIS, HEAT, THME, NINDPT, PATRAN 

COMMON /MICMAT/ EA(2) ,ET(2) ,GA(2) ,FNA(2) ,FNT(2) ,ALPA(2) ,ALPT(2) , 
& VPBP(6,2) ,VPR0B(2,12) 

COMMON /MICRO/ VF , MATNUM ( 2 , 2) , NST , NMT , NB , NG , IDP , NSFD , 

& NCMD 

COMMON /MICR02/ VF1 ,VC1 ,VF2 ,VC2 
COMMON /MICR03/ ICOUNT.NPLVL 
COMMON /MICR04/ A(NAM,NAM12) ,XH(7) ,XL(7) 

COMMON /MICR05/ CTEMP 

COMMON /MICR06/ SIGAL(181) ,DSIGAL(181) 


c /\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ 
c 


trVv.’ 


h- 


£k.. 


. IN FENTIO^ALLV 
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PRECEDING PAGE BLANK NOT FILMED 


DIMENSION D(NGENS.NGENS) ,G(NGENS) ,E(NGENS) ,DE(NGENS) ,S(NGENS) 
DIMENSION TEMP(l) ,DTEMP( 1) 

c \/\/\/\/\/\/\/\/\/\/\/\/\/\/\A/\A/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\ 
c 

C 

C DIMENSION STATEMENTS FOR GMC 

C 

C 

DIMENSION STRO(6) ,STRC(6) ,ASTRO(6) ,ASTRC(6) ,CD(6,6) ,CDI(6,6) 
DIMENSION STRAIN(6) .STRESS (6) 

DIMENSION EPS0G(6) ,EPIOG(6) ,STROG(6) 

DIMENSION EPS0SC(6) ,EPIOSC(6) ,STR0SC(6) 

DIMENSION EPSI(6) ,STRESR(6) ,DER(6) ,ALPHR(6) 

DIMENSION SA(26) ,PP(6) ,ZDRR(6) ,R(6) ,SD(6) ,EPIINC(6) 

DIMENSION IPVT(NAM) 

C 

DIMENSION VS(NSM) ,VS1(NSM,NMM) ,VSMR(NSM) ,MATT(NMM) 

C 

DIMENSION GAM(3 ,NMM) ,GAMS(3) ,ALPS(3) ,ALPSI(3) 

DIMENSION CS(6,6,NMM) ,CG(6,6) ,CI(6,6) ,CPDRV(6,6) 

DIMENSION SG(6) ,EG(6) ,EPG(6) ,ETG(6) ,H(6) ,U(6,NSM) 

DIMENSION SS(6) ,ES(6) ,EPS(NAM) ,ETS(NAM) ,UT(6,NSM) ,HT(3) 

DIMENSION SIGLC36) ,DSIGL(36) ,ES0(6) ,EPSO(6) ,ETSO(6) 

DIMENSION HILL(NSM,6,6) ,HILLT(NSM,3) 

C 

DIMENSION THRM(IO) 

c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA/ 

C 

c 

C Statement function: contracted notation indexation rule 
C for generalized method of cells 

c 

C 

PCON (IPCON , JPCON) = IPCON + JPCON 

MC0N(IMC0N, JMCON) = INT( 1/(ABS(IMC0N -JMCON)+l) ) 

NCON(INCON, JNCON) = MC0N(INC0N, JNCON)*PCON(INCON, JNC0N)/2 + 

& (1-MCON (INCON , JNCON) ) * (9-PCON (INCON , JNCON) ) 

c 

C GMC 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


Purpose: generalized cells model for metal matrix long-fiber 

composites the fibers are oriented in the xl-direction 


NBM = MAX. NUMBER OF SUBCELLS IN BETA-DIRECTION (X2) 

NGM = MAX. NUMBER OF SUBCELLS IN GAMMA-DIRECTION (X3) 

NMM = MAX. NUMBER OF MATERIALS 

NSM = NBM*NGM = MAX. NUMBER OF SUBCELLS IN A 

REPRESENTATIVE CELL 

CS = CONSTITUENT ELASTIC STIFFNESS TENSOR 

CG = GLOBAL (EFFECTIVE) ELASTIC STIFFNESS TENSOR 

Cl = INVERSE OF CG 

SG = GLOBAL STRESS TENSOR 

EG = GLOBAL STRAIN TENSOR 

EPG = GLOBAL PLASTIC STRAIN TENSOR 

ETG = GLOBAL THERMAL STRAIN TENSOR 

SS = SUBCELL STRESS TENSOR 

ES = SUBCELL STRAIN TENSOR 

EPS = SUBCELL PLASTIC STRAIN TENSOR 

EPD = SUBCELL PLASTIC STRAIN-RATE TENSOR 

ETS = SUBCELL THERMAL STRAIN TENSOR 

DTEMP = DEVIATION FROM THE REFERENCE TEMPERATURE 

NINDPT= FLAG WHEN MATERIAL CONSTANTS ARE INDEPENDENT 
OF TEMPERATURE 

NOTE THAT EG(4)=2*EPS23 ,EG(5)=2*EPS13 , EG(6)=2*EPS12 
AND SIMILAR RELATIONS HOLD FOR EPG(4) ,EPG(5) ,EPG(6) 


Limit is currently set for a maximum of 4 subcells ! ! ! ! 


I/O FILES: 

UNIT I CONTENTS 


#60 - GRID INFO 
#61 - PATGEO.GMC 
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C #62 - PATSTR . GMC 

C #66 - PROGRAM TRACE FILE 

C 

C 

C Set parameters for generalized method of cells 

C 

C 

NPLVL = 1 

C 

FORDIS = .TRUE. 

HEAT = .FALSE. 

THME = .FALSE. 

NINDPT = .TRUE. 

C 

IF(INC.EQ.O) THEN 
ICOUNT = 1 
ELSE 

ICOUNT = ICOUNT + 1 
END IF 
C 

IDP = 1 

NB =2 

NG =2 

PATRAN = .FALSE. 

NCMD = 1 
TEMPR = 0 . 

DTEMPR = 0 . 

C 

C Assign subcell material id. 

C 

MATNUM (1,1) = 1 
MATNUM(1 ,2) = 2 
MATNUM (2, 1) = 2 
MATNUM (2, 2) = 2 
MATT(l) = 1 
MATT (2) = 2 

C 

C 
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C Material constants: 

C NMT - number of different materials 
C NST - number of subcells 
C VF - fiber volume fraction 
C VC - coating volume fraction 

C 

C 

NMT = 2 
NST = 4 
VF = 0.35 
VC = 0. 

C 

C 

C Elastic constituent material properties: 
C 

C Fiber: 

C 

c 

EA(1) = 58. D3 

ET(1) = 58. D3 

GA(1) = 23.2D3 

FNA(l) = 0.25 
FNT(l) = 0.25 
ALPA(l) = 1. 

ALPT(l) = 1. 

C D_0 

VPR0B(1 , 1) = 0. 

C Z_0 

VPR0B(1 ,2) = 1. 

C Z_1 

VPR0B(1 ,3) = 1. 

C m 

VPR0B(1 ,4) = 1- 

C n 

VPR0B(1 ,5) = 1. 

C Q 

VPROBCl ,6) = 1. 

C 

C 
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C Matrix 
C 

c 

EA(2) = 11233. 

ET(2) = 11233. 

GA(2) = 4255. 

FNAC2) = 0.32 
FNT(2) =0.32 
ALPA(2) = 1. 

ALPTC2) = 1. 

C 

VPR0B(2,1) = 10. 

VPR0B(2,2) = 1600. 

VPR0B(2,3) = 1600. 

VPR0B(2,4) = 10. 

VPR0B(2,5) = 0.2 
VPR0B(2,6) = 10. 

C 

C 

C Form elastic material stiffness 

C 

C 

IF ( (NCYCLE+INC) .EQ. 0) THEN 
GOTO 5100 

ELSEIF (NCYCLE .EQ. 0) THEN 
GOTO 5100 
ELSE 
C 

END IF 

C— 

C Set global time step: GDT 

C 40 steps 

C eps.max = 0.02 rate = 0.001 

C - 

C 

GDT = 0.5 
C 

C 
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C -CTEMP(l)} - cycle number 
C {TEMP(2)> - global time step 

C 

C 

TEMP(l) = NCYCLE 
IF(INC .GT. 0) THEN 
IF (NCYCLE .EQ. 1) THEN 
DTEMPC2) = GDT 
END IF 
ELSE 

TEMP (2) = 0. 

ENDIF 

C 

c 

C {DE} - incremental strain vector 
C 

C GMC global strain rate vector {EG} 

C — 

C 

LIM = 37 + NST*36 
DO 1002 K = l.LIM 
SIGAL(K) = TEMP(K+2) 

DSIGAL(K) = DTEMPCK+2) 

1002 CONTINUE 
C 

DO 1003 J=1 ,6 

EG( J) = DE( J)/GDT 
DSIGAL(J) = EG(J) 

EPSOG(J) = TEMPCJ+2) 

STROG(J) = TEMP( J+8) 

EPIOG(J) = TEMPCJ+14) 

STRESS (J)= S(J) 

1003 CONTINUE 
C 

C Begin generalized method of cells code 
C 

C \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ 
C 
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5100 CONTINUE 
C 

C 

C Square packing (original method of cells with 4 subcells) 

C 

C 

IFCIDP.EQ.1 ) THEN 
XH ( 1 ) =SQRT ( VF) 

XH(2)=1-SQRT(VF) 

XL(1)=XH(1) 

XL(2)=XH(2) 

C 

C — 

C Triangle packing 

C 

C 

ELSEIF(IDP .EQ. 2) THEN 
IFCVF.LT. 0.288675) THEN 

XA=SQRT(SQRT(3 . DO) *VF/2 . DO) 

XB=(l-XA)/2.0 
XC=SQRT(3 . DO/4 . DO) -XA 
XH(l)=XB-XA/2 
XH(2)=XA 
XH(3)=XB-XA/2 
XH(4)=XA 
XL(1)=XA 
XL(2)=XC 
XL(3)=XA 
XL(4)=XC 
ELSE 

XA=DSQRT(DSQRT(3 . DO) *VF/2 . DO) 

XB=(l-XA)/2 

XC=DSQRT (3 . DO/4 . DO) -XA 

XH(1)=2.0*XB 

XH(2)=XA/2 .O-XB 

XH(3)=2.0*XB 

XH(4)=XA/2.0-XB 

XL(1)=XA 

XL(2)=XC 


34 


XL(3)=XA 
XL(4)=XC 
END IF 
C 

C 

C Square diagonal packing 

C 

C 

ELSEIFCIDP .EQ. 3) THEN 
XA=DSQRT (VF/2 . DO) 

XH(l)=(l.-2.*XA)/2. 

XH(2)=XA 

XH(3)=(l-2*XA)/2 

XH(4)=XA 

XL(1)=XA 

XL(2)=(l-2*XA)/2 

XL(3)=XA 

XL(4)=(l-2*XA)/2 

C 

C — 

C Cross shape fiber (square packing) 

C 

C 

ELSEIFCIDP .EQ. 4) THEN 

XB=-2 . *XA+DSQRT(4 . D0*XA**2 . D0*VF) 

XH(1)=XA 
XH(2)=XB 
X H(3)=XA 
XH(4)=1-2*XA-XB 
XL(1)=XH(1) 

XL(2)=XH(2) 

XL(3)=XH(3) 

XL(4)=XH(4) 

C 

C 

C Coated fiber (square packing); vc=coating volume fraction 

C 

C 

ELSEIFCIDP .EQ. 5) THEN 
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XA=DSQRT(VF) 

XH(2)=XA 

XH ( 1 ) = ( -XA+DSQRT (VF+VC) ) / 2 . 

XH(3)=XH(1) 

XH(4)=1-(XH(1)+XH(2)+XH(3) ) 

XL(l)=XH(l) 

XL(2)=XH(2) 

XL(3)=XH(3) 

XL(4)=XH(4) 

C 

c 

C Circular fiber approximation 
C (49 cell model) 

C 

C 

ELSEIFCIDP .EQ. 6) THEN 
RADIUS=DSQRT ( VF/PAI ) 

HH=1 .0 

XH (2) =DSQRT (PAI) *RADIUS/DSQRT (52 . DO) 
XH(3)=XH(2) 

XH(4)=4.0*XH(2) 

XH(5)=XH(3) 

XH(6)=XH(2) 

XH(l)=(HH-2.0*XH(2)-2.0*XH(3)-XH(4))/2.0 

XH(7)=XH(l) 

XL(l)=XH(l) 

XL(2)=XH(2) 

XL(3)=XH(3) 

XL(4)=XH(4) 

XL(5)=XH(5) 

XL(6)=XH(6) 

XL(7)=XH(7) 

END IF 
C 

c 

C Calculate subcell volume 
C {VS} 

C 

C 
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SUMH=0 

DO 1010 1=1, NB 
SUMH=SUMH+XH ( I ) 

1010 CONTINUE 
SUML=0 

DO 1012 1=1, NG 
SUML=SUML+XL(I) 

1012 CONTINUE 
C 

DO 1113 IB=1 ,NB 
DO 1114 IG=1,NG 

VS(NG*(IB-1)+IG) = XH(IB)*XL(IG) 

1114 CONTINUE 

1113 CONTINUE 
C 

C 

C Subcell information output 
C (optional) 

C 

C 

IF (NPLVL . GE . 5) THEN 
WRITE (66, 683) 

WRITE (66, 684) 

END IF 
C 

TOTV=0.0 
DO 1120 IB=1,NB 
DO 1120 IG=1 ,NG 
IS=NG*(IB-1)+IG 
MN =M ATNUM ( IB , I G ) 

KK = MATT(MN) 

TQTV=TQTV + VS (IS) 

IF (NPLVL. GE.S) THEN 

WRITE (66 , 685) IB , IG , IS , MN , KK , VS (IS) 
ENDIF 

1120 CONTINUE 
C 

c 

C Calculate subcell material volume ratio 
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C -CVSMR} 

C 

C 

DO 1130 MN=1 ,NMT 
VSMR(MN)=0 
DO 1130 IS=1,NB*NG 
VS1(IS,MN)=0 
1130 CONTINUE 

DO 1140 IB=1 ,NB 
DO 1140 IG=1 ,NG 
IS=NG*(IB-1)+IG 
MN=MATNUM(IB,IG) 

VS1 (IS,MN)=VS(IS) 

1140 CONTINUE 

DO 1141 MN=1 ,NMT 
DO 1141 IS=1 ,NB*NG 

VSMR (MN ) =VSMR (MN ) +VS 1 ( IS , MN ) 

1141 CONTINUE 

DO 1142 MN=1 ,NMT 

VSMR (MN) =VSMR (MN ) /TOTV 

IF(NPLVL.GE.5) WRITE(66 , 143) MN.VSMR(MN) 

1142 CONTINUE 
C 

C 

C Initialization 

c 

c 

NST=NB*NG 
DO 1150 1=1,6 
DO 1150 >1,6 
DO 1150 K=1 ,NMM 
CS(I , J,K)=0 
1150 CONTINUE 

DO 1160 1=1,6 
DO 1160 J=1 ,6 
CI(I,J)=0 
CG(I , J)=0 
1160 CONTINUE 
C 


38 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DEFINE THE TEMPERATURE-DEPENDENT PROPERTIES AT 6 DIFFERENT 
TEMPERATURES: THRM(I) ,1=1... 6 

FOR INTERMEDIATE TEMPERATURES APPLY A LINEAR INTERPOLATION 
===> CURRENTLY ONLY FOR BODNER MODEL <=== 


IF (HEAT .OR. THME) THEN 
IF(NCMD.EQ.l) THEN 
TEMPR = CTEMP 

IF(TEMPR.GE.THRM(1) .AND.TEMPR.LE.THRM(2)) THEN 
1=2 

ELSEIF (TEMPR . GE . THRM (2) . AND . TEMPR . LE . THRM(3) ) THEN 
1=3 

ELSEIF (TEMPR. GE. THRM (3) . AND. TEMPR. LE. THRM (4)) THEN 
1=4 

ELSEIF (TEMPR . GE . THRM (4) . AND . TEMPR . LE . THRM(5) ) THEN 
1=5 

ELSEIF (TEMPR . GE . THRM ( 5 ) . AND . TEMPR . LE . THRM ( 6) ) THEN 
1=6 


& 


& 

& 

& 


ELSE 

WRITE (66,*) ’GMC: **** ERROR ****’ 

WRITE(66,*) * TEMPERATURE IS OUT OF RANGE ’, TEMPR 
STOP 
END IF 

DIFF= (TEMPR-THRM (I))/ (THRM ( I ) -THRM ( I - 1 ) ) 

DO 1200 KK=1 ,NMT 

IF(NPLVL.GE.5) WRITE(66,*) ’GMC: NM ’ ,NM, ’ KK ’ ,KK 

EA(KK) = EATM(I ,KK)+ 

(EATM(I,KK)-EATM(I-1,KK))*DIFF 
FNA(KK) = FNATM(I ,KK)+ 

(FNATM(I ,KK)-FNATM(I-1 ,KK) )*DIFF 
ET(KK) = ETTM(I ,KK)+ 

(ETTM (I , KK) -ETTM (I- 1 , KK) ) *DIFF 
FNT(KK) = FNTTM(I ,KK)+ 

(FNTTM(I ,KK)-FNTTM(I-1 ,KK))*DIFF 
GA(KK) = GATM(I ,KK)+ 

(GATM(I ,KK)-GATM(I-1 ,KK) )*DIFF 
ALPA(KK) = ALPATM(I ,KK)+ 
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c 

c 

c 

C1200 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cl210 


(ALPATM(I ,KK)-ALPATM(I-1 ,KK) )*DIFF 
ALPT(KK) = ALPTTM(I ,KK)+ 

(ALPTTM(I ,KK)-ALPTTH(I-1 ,KK) )*DIFF 

CONTINUE 

ELSEIF (NCMD . NE . 1) THEN 
CALL CONEVL ( 

DO 1210 KK=1 ,NMT 


EA(KK) 

ET(KK) 

GA(KK) 

FNA(KK) 

FNT(KK) 

ALPA(KK) 

ALPT(KK) 

CONTINUE 

ELSE 

WHITE (66,*) 

STOP 
END IF 
END IF 


PEM(l.KK) 

PEM(2,KK) 

PEH(3,KK) 

PEM(4,KK) 

PEM(5,KK) 

PEM(6,KK) 

PEM(7,KK) 


’GMC: ERROR IN TEMPERATURE ’ 
’ DEPENDENT CONSTANTS’ 


C 

C- 

c 

c- 

c 


Elastic stiffness matrix [CS] for each subcell 


DO 1230 NM= 1 , NMT 

GT = 0.5*ET(NM)/(1+FNT(NM)) 

FK = 0.25*EA(NM)/(0.5*(1-FNT(NM) )* 

(EA (NM) /ET(NM) ) -FNA(NM) **2) 
CS(l.l.NM) = EA(NM) + 4*FK*FNA(NM)**2 
CS(2 , 1 ,NM) = 2*FK*FNA(NM) 

CS(3,1,NM) = 2*FK*FNA(NM) 

CS(1 ,2 ,NM) = 2*FK*FNA(NM) 

CS(2,2,NM) = FK+GT 
CS(3,2,NM) = FK-GT 
CS(1 ,3 ,NM) = 2*FK*FNA(NM) 

CS(2,3,NM) = FK-GT 
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CS(3,3,NM) = FK+GT 
CS(4,4,NM) = GT 
CS(5,5,NM) = GA(NM) 

CS(6,6,NM) = GA(NM) 

C 

IFCNPLVL.GE.5) THEN 

c IF(NPLVL.GE.2 .and. inc.ge.O .and. nn.eq.l) then 

write(66,*) ’nm ’ ,nm 
write(66,*) ’ns ’ ,ns 

WRITE (66, 8266) 

WRITE (66, 8267) CS(1,1,NM) ,CS(1,2,NM) ,CS(1,3,NM) 
WRITE (66 , 8267 ) CS(2,1,NM) ,CS(2,2,NM) ,CS(2,3,NM) 
WRITE (66 , 8267 ) CS(3,1,NM) ,CS(3,2,NM) ,CS(3,3,NM) 
WRITE (66, 8268) CS(4,4,NM) ,CS(5,5,NM) ,CS(6,6,NM) 

END IF 
C 

C 

C Thermal stress vectors [GAM] of subcells 

C 

C 

GAM(1 ,NM)=CS(1 , 1 ,NM)*ALPA(NM)+(CS(1 ,2 ,NM)+ 

& CS(1 ,3,NM))*ALPT(NM) 

GAM(2,NM)=CS(1 ,2,NM)*ALPA(NM)+(CS(2,2,NM)+ 

& CS(2,3,NM))*ALPT(NM) 

GAM (3 , NM) =CS (1,3, NM) *ALPA (NM) + (CS (2 , 3 ,NM) + 

& CS (3 , 3 , NM) ) *ALPT(NM) 

1230 CONTINUE 
C 

C 

C Compute the concentration factor 
C initialization of the a-matrix 

C - 

C 

DO 1240 J=l, 6*NST + 6 + 6*NST 
DO 1250 1=1 ,6*NST 
A(I , J)=0 . 000 
1250 CONTINUE 

1240 CONTINUE 
C 
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NSR = NG*(IB-1) + IGR 
MSL = MATNUM(IB.IGL) 

MSR = MATNUM(IB.IGR) 

NA=NA + 1 
DO 1650 NC=1 ,6 

A(NA,6*(NSL-1)+NC) = CS(MC,NC,MSL) 
A(NA,6*(NSR-1)+NC) = -CS(MC,NC,MSR) 
A(NA ,6*(NST+NSL)+NC) = CS(MC ,NC ,MSL) 
A(NA,6*(NST+NSR)+NC) = -CS(MC ,NC ,MSR) 
1650 CONTINUE 

1660 CONTINUE 

1670 CONTINUE 

1680 CONTINUE 
C 

C 

C Displacement continuity conditions 

C 

C 

C ES11(B,G)=EG11 B=1 .... ,NB ; G=1,...,NG 

C 

C 

DO 1710 IB=1 ,NB 
DO 1705 IG=1 ,NG 
NA=NA+1 

A(NA ,6*NB*NG+1) = 1 
NS = NG*(IB-1)+IG 
A(NA ,6*(NS-1)+1) = 1 
1705 CONTINUE 

1710 CONTINUE 
C 

C 

C SUM [H(B)*ES22(B,G)]=H*EG22 G=1,...,NG 

C 

C 

DO 1720 IG=1 ,NG 
NA=NA+1 

A(NA ,6*NB*NG+2) = SUMH 
DO 1715 IB=1 ,NB 
NS = NG*(IB-1)+IG 
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A(NA ,6*(NS-l)+2) = XH(IB) 


1715 

CONTINUE 


C 



C 

C SUM 

[H(B)*ES12(B,G)]=H*EG12 G=l, . . 

. ,NG 

c 



c 

NA=NA+1 

A(NA,6*NB*NG+6) = SUMH 
DO 1718 IB=1,NB 



NS = NG*(IB-1)+IG 
A(NA,6*(NS-l)+6) = XH(IB) 


1718 

CONTINUE 


1720 

CONTINUE 


C 



C 

C SUM 

[L(G)*ES33(B,G)]=L*EG33 B=l... 

. ,NB 

c 



c 

DO 1730 IB=1 ,NB 



NA=NA+1 

A(NA,6*NB*NG+3) = SUML 
DO 1725 IG=1 ,NG 



NS = NG*(IB-1)+IG 
A(NA,6*(NS-l)+3) = XL(IG) 


1725 

CONTINUE 


C 



c 

C SUM 

[L(G)*ES13(B,G)]=L*EG13 B=l,.. 

. ,NB 


c 

NA=NA+1 

A(NA,6*NB*NG+5) = SUML 
DO 1728 IG=1 ,NG 
NS = NG*(IB-1)+IG 
A(NA,6*(NS-l)+5) = XL(IG) 
1728 CONTINUE 

1730 CONTINUE 
C 
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c 

C SUM SUM [H (B) *L (G) *ES23 (B , G) ] =H*L*EG23 

C 

C 

NA=NA+1 

DO 1740 IB=1 ,NB 
DO 1735 IG=1 ,NG 

A(NA,6*NB*NG+4) = SUMH*SUML 
NS = NG*(IB-1)+IG 
A (N A, 6* (NS- 1) +4) = XH(IB)*XL(IG) 

1735 CONTINUE 

1740 CONTINUE 
C 

CALL SGEFA ( A , NAM , 6*NB*NG , IPVT , INFO) 

IF(INFO.NE.O) STOP 

DO 2110 1=1 ,6+6*NB*NG 

CALL SGESL ( A , N AM , 6*NB*NG , IPVT , A ( 1 , 6*NB*NG+I ) ,0) 

2110 CONTINUE 
C 

C 

C Compute the effective elastic stiffness matrix [CG] 

C 

C 

DO 2800 1=1,6 

DO 2790 J=1 ,6 
DO 2780 IB=1 ,NB 
DO 2770 IG=1 ,NG 
IS=NG*(IB-1)+IG 
IM=MATNUM(IB,IG) 

TMP=0 . 000 
DO 2760 M=1 ,6 

TMP=TMP + CS(I ,M ,IM)*A(6*(IS-1)+M,6*NST+J) 
2760 CONTINUE 

CG(I , J)=CG(I , J) + VS(IS)*TMP/(SUMH*SUML) 

2770 CONTINUE 

2780 CONTINUE 

2790 CONTINUE 
2800 CONTINUE 
C 
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c 

C MARC control statements 

C— 

IF ( (NCYCLE+INC) .EQ. 0) THEN 
CALL PLACE(CG,D,36) 

TEMP(l) = 0. 

TEMP (2) = 0. 

CALL ZEROR(SIGAL, 181) 

CALL ZEROR(DSIGAL , 181) 

CALL ZEROR(SIGL ,36) 

CALL ZER0R(DSIGL,36) 

RETURN 

ELSEIFCNCYCLE .EQ. 0) THEN 
CALL PLACE(CG,D ,36) 

RETURN 

ELSE 

CALL PLACE(CG,D,36) 

CALL PLACE(CG,CD,36) 

END IF 
C 

C 

C Invert the effective elastic stiffness matrix [CG] 

C 

C 

DET=CG(1 , 1)*(CG(2,2)*CG(3 ,3)-CG(2 ,3)**2) 

& -CG(1 ,2)*(CG(1 ,2)*CG(3,3)-CG(2,3)*CG(1 ,3) ) 

& +CG(1 ,3)*(CG(1 ,2)*CG(2,3)-CG(2,2) *CG(1 ,3) ) 

CI(1 , 1)=(CG(2,2)*CG(3 ,3)-CG(2,3)**2)/DET 
Cl (1 ,2)=(-CG(l ,2)*CG(3 ,3)+CG(l ,3)*CG(2,3))/DET 
CI(1 ,3)=(CG(1 ,2)*CG(2 ,3)-CG(l ,3)*CG(2,2) )/DET 
CI(2,2)=(CG(1,1)*CG(3,3)-CG(1,3)*CG(1,3))/DET 
CI(2,3)=(-CG(1 ,1)*CG(2,3)+CG(1 ,2)*CG(1 ,3))/DET 
Cl (3 ,3)=(CG(1 , 1)*CG(2 ,2) -CG(1 ,2) *CG(1 ,2) )/DET 
CI(4,4)=1/CG(4 ,4) 

CI(5,5)=1/CG(5,5) 

CI(6,6)=1/CG(6,6) 

CI(2,1)=CI(1,2) 

CI(3, 1)=CI(1 ,3) 

CI(3,2)=CI(2,3) 
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c 

IF(NPLVL.GE.5) THEN 
WRITE (66, 869) 

WRITE(66 ,91) 1. /Cl (1,1) 

WRITE (66 ,92) -CI(1,2)/CI(1,1) 

WRITE (66 ,93) l./CI(2,2) 

WRITE (66 ,94) -CI(2,3)/CI(2,2) 

WRITE(66 ,95) l./CI(3,3) 

WRITE(66 ,96) l./CI(4,4) 

WRITE (66 ,97) l./CI(5,5) 

WRITE(66 ,98) l./CI(6,6) 

ENDIF 

C 

C 

C Compute the effective coefficients 

C 

C 

DO 2801 1=1,3 
GAMS(I)=0. 

DO 2801 IG=1 ,NG 
DO 2801 IB=1 ,NB 
IS=NG*(IB-1)+IG 
MN=MATNUM(IB,IG) 

SUM=0 

DO 2761 K-1,3 

SUM=SUM +A(6*(IS-1)+K,6*NST+I)*GAM(K,MN) 

2761 CONTINUE 

GAMS (I ) =GAMS (I) +VS ( IS) *SUM/ (SUMH*SUML) 

2801 CONTINUE 
C 

ALPS(1)=CI(1,1)*GAMS(1)+CI(1,2)*GAMS(2)+CI(1,3)*GAMS(3) 
ALPS(2)=CI(2, 1)*GAMS(1)+CI(2,2)*GAMS(2)+CI(2,3)*GAMS(3) 
ALPS(3)=CI(3, 1)*GAMS(1)+CI(3,2)*GAMS(2)+CI(3,3)*GAMS(3) 
C 

C 

C Determine the effective cte directly from the model 

C 

C 

IF (HEAT .OR. THME) THEN 
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DO 3130 IB=1 ,NB 
DO 3130 IG=1 ,NG 
NS=NG*(IB-1)+IG 
NM=MATNUM(IB,IG) 

KK = NM 

ETS (6* (NS- 1) +1 ) =ALPA (KK) 

ETS (6* (NS- 1 ) +2) =ALPT (KK) 

ETS (6* (NS- 1 ) +3) =ALPT (KK) 

ETS (6* (NS- 1) +4) =0 
ETS (6* (NS- 1) +5) =0 
ETS (6* (NS-l)+6)=0 
C 

DO 3140 JJ=1 ,6 

DSIGAL(37+(36*(NS-1) )+( JJ+30))=ETS(6*(NS-1)+JJ)*DTEMPR 
3140 CONTINUE 

3130 CONTINUE 
C 

C 

c [UT] = [AP (B , G) ] *{ETS} 

C 

C 

DO 3145 NS=1 ,NST 
DO 3145 1=1,6 
UT(I ,NS)=0 
DO 3145 J=1 ,6*NST 

UT(I ,NS)=UT(I ,NS)+A(6*(NS-1)+I ,6*NST+6+J)*ETS( J) 

3145 CONTINUE 
C 

C 

C {HT}=SUM SUM (VS(B,G)*[CS]*([UT(B,G)]-[ETS(B,G)]))/TOTV 

C 

C 

DO 3150 1=1,3 
HT(I)=0 

DO 3150 IB=1 ,NB 
DO 3150 IG=1 ,NG 
NS=NG*(IB-1)+IG 
NM=MATNUM ( IB , I G ) 

DO 3150 M=1 ,3 
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HT(I)=HT(I)+VS(NS)* 

& CS (I , H , NH) * (UT(M , NS) -ETS (6* (NS-1) +M) ) / (SUMH*SUML) 

3150 CONTINUE 
C 

C 

C {ALPS}=- [Cl] *{HT> 

C 

c 

DO 3160 1=1,3 
ALPS (I) = 0. 

DO 3160 M=1 ,3 

ALPS (I) = ALPS (I) -Cl (I ,M)*HT(M) 

3160 CONTINUE 
C 

C 

C CALCULATE INSTANTANEOUS THERMAL EXPANSION 

C 

C 

c CALL ALPEVL ( ALPS , VS , SUMH , SUML , ALPSI ) 

ELSE 

CALL ZEROR(ETS.NAM) 

ENDIF 

C 

C 

C Determination of hill concentration 
C tensors [HILL] in subcell (B,G) 

C [HILL (B , G) ] = [CS (B , G) ] * [A (B , G) ] * [Cl] 

C 

C 

DO 3170 IB=1 ,NB 
DO 3170 IG=1 ,NG 
NS=NG*(IB-1)+IG 
NM=MATNUM(IB,IG) 

DO 3170 1=1,6 
DO 3170 J=1 ,6 
HILL(NS,I, J)=0 
DO 3170 M=1 ,6 
DO 3170 N=1 ,6 

HILL (NS , I , J)=HILL (NS , I , J) + 
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CS (I , M , NM) *A (6* (NS- 1) +M , 6*NST+N) *CI (N , J) 


& 

3170 CONTINUE 
C 

IFClNC.ge.O .AND. NPLVL.GE.5) THEN 
WRITE (6, 7 129) 

DO 3180 IB=1 ,NB 
DO 3180 IG=1 ,NG 
NS=NG*(IB-1)+IG 
WRITE(6,7126) IB.IG 
DO 3181 1=1,6 

WRITE(66,7128) (HILL (NS, I ,J) , J=1 ,6) 

3181 CONTINUE 
3180 CONTINUE 
END IF 
C 



C Determination of hill thermal concentration 
C vector [HILLT] in the subcells 

C [HILLT(B,G)]=-[CS(B,G)]*([UT(B,G)]-[ETS(B,G)] + [A(B,G)]*{ALPS» 

C 

C 

IF (HEAT .OR. THME) THEN 
DO 3185 1=1,3 
DO 3185 IB=1 ,NB 
DO 3185 IG=1 ,NG 
NS=NG*(IB-1)+IG 
NM=MATNUM ( IB , I G ) 

HILLT(NS,I)=0 
DO 3185 M=1 ,3 
HT(H)=0 
DO 3190 N=1 ,3 

HT(M)=HT(M)+A(6*(NS-1)+M,6*NST+N)*ALPS(N) 

3190 CONTINUE 

HILLT (NS , I ) =HILLT (NS , I ) 

& -CS (I , M , NM) * (UT (M , NS) -ETS (6* (NS- 1 ) +M) +HT (M) ) 

3185 CONTINUE 
C 

IF(INC.EQ.O .AND. NPLVL.GE.5) THEN 
WRITE (66, 7139) 
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DO 3200 IB=1 ,NB 
DO 3200 IG=1 ,NG 
NS=NG*(IB-1)+IG 

WRITE (66 ,7138) IB ,IG, (HILLT(NS ,1) ,1=1 ,3) 

3200 CONTINUE 

ENDIF 
ENDIF 
C 

C 

C Finish elastic calculations 

C 

C Begin inelastic calculations 

C 

C 

DO 8800 IB=1 ,NB 
DO 8750 IG=1 ,NG 

NS = NG*(IB-1) + IG 
NM = MATNUM(IB ,IG) 

C 

DO 2222 J=1 ,6 

EPSOSC(J) = TEMP ((2 + 37 + (36*(NS-1))) +J) 

STROSC(J) = TEMP ( (2 + 37 + (36*(NS-1))) +J+6) 

EPIOSC(J) = TEMP ( (2 + 37 + (36*(NS-1))) +J+12) 

2222 CONTINUE 

C 

C 

C Copy subcell quantities from 
C {SIGAL} TO -CSIGL} AND 
C -CDSIGAL} TO {DSIGL} 

C 

C 

DO 7705 KK=1 ,36 

SIGL(KK) = SIGAL(37+(36*(NS-1))+KK) 

DSIGL (KK) = DSIGAL (37 + (36* (NS- 1) ) +KK) 

7705 CONTINUE 

C 

C 

C Call constitutive model to get: 

C inelastic strain rate {EPS} 
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C and state variable rates 

C 

C 

IF(NCMD .EQ. 1) THEN 

CALL BODNER(DSIGL , SIGL , STROSC , NM , NS , INC , N , NN) 
ELSE 
C 

C Insert other viscoplastic material models here 
C 

WRITE(66,*) ’ ILLEGAL MODEL’ 

STOP 
END IF 
C 

C 

C Copy subcell quantities into appropriate 
C postions in {SIGAL} and {DSIGAL} 



C 

DO 7700 KK=1 ,36 

SIGAL(37+(36*(NS-1) )+KK) = SIGL(KK) 
DSIGAL(37+(36*(NS-1) )+KK) = DSIGL(KK) 

7700 CONTINUE 

C 

C 

C Copy subcell plastic strain rate {EPS} 

C from {DSIGAL} 



C 

EPS(6*(NS-1)+1)=DSIGAL(37+(36*(NS-1)+13) ) 

EPS (6* (NS-1 ) +2) =DSIGAL (37+ (36* (NS- 1) +14) ) 

EPS (6* (NS- 1 ) +3) =DSIGAL (37+ (36* (NS- 1 ) + 15) ) 
EPS(6*(NS-1)+4)=DSIGAL(37+(36*(NS-1)+16)) 

EPS (6* (NS-1) +5) =DSIGAL (37+ (36* (NS- 1 ) +17) ) 

EPS (6* (NS-1 ) +6) =DSIGAL (37+ (36* (NS- 1 ) + 18) ) 

8750 CONTINUE 
8800 CONTINUE 

C 

C End subcell loop 
C 
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c 

c 

C Compute the average composite plastic strain rate {EPG} 

C 

C 

DO 8930 IS=1 ,NST 
DO 8920 1=1,6 
U(I,IS)=0.000 
DO 8910 J=1,6*NST 

U(I,IS)=U(I,IS) + A(6*(IS-l)+I,6*NST+6+J)*EPS(J) 

8910 CONTINUE 

8920 CONTINUE 

8930 CONTINUE 
C 

DO 8950 1=1,6 
H(I)=O.ODO 
DO 8960 IB=1 ,NB 
DO 8960 IG=1 ,NG 
IS=NG*(IB-1)+IG 
IM=MATNUM(IB ,IG) 

DO 8940 J=1 ,6 

H(I)=H(I) +VS(IS)*CS(I , J ,IM)*( U(J ,IS)-EPS(6*(IS-1)+J) ) 
& / (SUMH*SUML) 

8940 CONTINUE 

8960 CONTINUE 

8950 CONTINUE 

C 

DO 8980 1=1,6 
EPG(I)=O.ODO 
DO 8970 >1,6 

EPG(I) = EPG(I) - Cl (I , J)*H(J) 

8970 CONTINUE 

DSIGALCI+12) = EPG (I) 

8980 CONTINUE 

C 

C 

C Global strains {DSIGAL(l-6)} 

C 

c 
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C Subcell loop 

c 

C 

DO 9775 IB=1 ,NB 
DO 9775 IG=1 ,NG 
C 



C Compute the average strain rate {ES> in the subcell 
C {ES}= [A(B,G)]*{EG> + [AP(B ,G)] *{EPS> + [AP(B ,G)] *{ETS> 

C - 

C 

NS = NG*(IB-1) + IG 
NM = MATNUM(IB,IG) 

DO 9150 M=1 ,6 

ES(M) = U(M,NS) + UT(M,NS)*DTEMPR 
DO 9100 N=1 ,6 

ES(M) = ES(M) + A(6*(NS-1)+M,6*NST+N)*EG(N) 

9100 CONTINUE 

DSIGAL(37+(36*(NS-1))+M) = ES(M) 

9150 CONTINUE 

C 

C 

C Compute the average stress rate {SS} in the subcell 
C {SS}= [CS] * (ES (B , G) -EPS (B , G) -ETS (B , G) ) 

C 

C 

CALL ZER0R(SS ,6) 

DO 9250 M2=l ,6 
DO 9200 N2=l ,6 

SS(M2) = SS(M2) + CS(M2,N2,NM)*CES(N2)-EPS(6*(NS-1)+N2)- 
& ETS(6*(NS-1)+N2)*DTEMPR) 

9200 CONTINUE 

C 
C 
C 

IF(NINDPT) GOTO 9010 
C 

C 

C Add addtional term in stress rate equation 
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C when [CS] is a function of temperature 

C 

C 

IF (HEAT .OR. THME) THEN 
IFCM2.EQ.1) THEN 

c CALL CSEVAL (CPDRV) 

ENDIF 

IFCM2.EQ.1) THEN 
DO 9701 1=1,6 

ESO(I) = SIGAL(37+(36*(NS-1) )+I) 

EPSO(I) = SIGAL(37+(36*(NS-1)+(I+12))) 

ETSO(I) = SIGAL(37+(36*(NS-l)+(I+30) ) ) 

TEMPO = SIGAL(37) 

9701 CONTINUE 

ENDIF 

DO 9251 N3=l ,6 

ADDTRM = CPDRV (M ,N3) *DTEMPR* (ESO (N3) -EPS0(N3) -ETSO (N3) ) 
SS(M2) = SS(M2) + 

& CPDRV (M2 , N3)*DTEMPR* (ESO (N3) -EPSO (N3) -ETSO (N3) ) 

9251 CONTINUE 

ENDIF 
C 

9010 CONTINUE 

C 

C 

C Copy subcell stress rate {SS> to {DSIGAL} 

C 

C 

DSIGAL(37+(36*(NS-1) )+(M2+6) ) = SS(M2) 

9250 CONTINUE 
9775 CONTINUE 
C 

C — — 

C Compute the average composite 
C thermal strain rates -CETG} 

C {ALPSI} - instantaneous thermal 
C expansion coefficient 

C - 

C 
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CALL ZER0R(ETG,6) 

IF (HEAT .OR. THME) THEN 
ETG(l) = ALPSI ( 1 ) *DTEMPR 
ETG(2) = ALPSI (2) *DTEMPR 
ETG(3) = ALPSI (3) *DTEMPR 
ETG(4) = 0 
ETG(5) = 0 
ETG(6) = 0 
END IF 
C 

C 

C Copy global thermal strain rates 
C to {DSIGAL (31-36)} 

C 

C 

DO 9889 JJ=1,6 

c dsigal(j j+30) = etg(jj) 

DSIGAL (JJ+30) =0.0 

9889 CONTINUE 
C 

C 

C Compute the average composite stress rate field {SG} 

C and copy to {DSIGAL(7-12)} 

C 

C 

DO 9890 1=1,6 
SG(I)= 0.000 
DO 9885 J=1 ,6 

SG(I) = SG(I) + CG(I , J)*(EG( J)-EPG( J)-ETG( J) ) 

9885 CONTINUE 

DSIGAL(I+6) = SG(I) 

9890 CONTINUE 
C 

c /\ /\ /\ A A A A A A A A A A A A A A A A A A A A 
c 

C End generalized method of cells code 
C 

C 
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C Integrate rates using forward euler 

C — — 

C 

C 

C Form MARC increment vector {DTEMP} 

C — 

C 

LIM = 37 + NST*36 
DO 5000 K = l.LIM 

DTEMPCK+2) = GDT*DSIGAL(K) 

5000 CONTINUE 
C 

C 

C {G> - MARC vector of inelastic strain increment 

C — 

C 

CALL ZER0R(G ,6) 

DO 5010 I = 1,6 
DO 5010 J = 1,6 

G(I) = G(I) + GDT*CD(I , J)*DSIGAL(J+12) 

5010 CONTINUE 
C 

C 

C Update MARC stress vector {S} 

C from {SIGAL(7-12» 

c 

c 

CALL ZER0R(S ,6) 

DO 5020 I = 1,6 

S(I) = SIGALU+6) + GDT*DSIGAL(I+6) 

5020 CONTINUE 
C 

RETURN 

C 

C Format statements 
C 

91 F0RMATC5X, ’Ells=’ ,E11.3) 
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92 F0RMAT(5X,’N12s=\E11.3) 

93 F0RMAT(5X, ’E22S=’ ,E11.3) 

94 FORMAT (5X,’N23S=’ ,E11.3) 

95 FORMAT(SX,’E33S=’,E11.3) 

96 FORMAT(5X, ’G23S=’ , El 1.3) 

97 FORMAT(5X, ’G13S=’ ,E11.3) 

98 F0RMAT(5X, ’G12S=’ ,E11 .3) 

99 FORMAT(5X, ’ALPHA11S=’ ,E11 .3) 

C 

143 F0RMAT(10X, ’Material No. =’ ,12, 7X, ’Volume Ratio = ’.E10.3,/) 

C 

683 F0RMAT(10X, ’Subcell Data:’//) 

684 FORMAT ( IX , ’ (BETA I GAMMA) ’ ,2X, ’SC #’ ,2X, ’SC MAT. #’,2X, 
ft ’SC MAT. TYPE’ ,2x, ’SUBCELL VOLUME’/) 

685 FORMAT (2X , 12 , IX , ’ I ’ ,I2,5X,I3,5X,I2,5X,I2,12X,E10.3) 

C 

866 F0RMAT(///5X,’CG EFFECTIVE STIFFNESS MATRIX’//) 

867 F0RMAT(5X,E10.3,3X,E10.3,3X,E10.3) 

868 FORMAT (44X ,E10 . 3/ 

ft 57X.E10.3/ 

ft 70X.E10.3) 

869 F0RMATC//5X, ’EFFECTIVE ENGINEERING MODULI’//) 

870 FORMAT (5X, ’EFFECTIVE THERMAL EXPANSION COEFFICIENTS’//) 

C 

990 F0RMAT(5X, ’ ALPHA22S=’ ,E11 .3) 

991 FORMAT C5X,’ALPHA33S=’ ,E11.3) 

C 

6333 FORMAT (5X, ’EG’ ,1X,6E10.3,/) 

6334 FORMAT (5X, ’SG’ ,1X,6E10.3,/) 

6335 FORMAT (5X,’EPG’ .1X.6E10.3,/) 

6336 F0RMAT(5X , ’ ETG ’ , IX , 6E10 . 3 , /) 

C 

7134 FORMAT (/,5X, ’The Effective Coefficients Of Thermal Expansion’,/, 

& 5X, ’Obtained Directly From The Model (And Not Via Levin Formula)’ 
ft /,10X,3E11.3,/) 

7126 FORMATC/ , 10X , ’ IB=’ ,I3,2X, ’IG=’ ,13) 

7128 FORMAT (5X,6E1 1.4) 

7129 F0RMAT(19X, ’Hill Concentration Tensors Of The Subcells’) 

7138 FORMAT (5X, ’IB=’ ,12, 2X, ’IG=’ ,I2,2X,3E11 .4) 
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7139 F0RMAT(/,7X, ’Hill Subcell Thermal Concentration Vectors’,/) 
C 

8266 F0RMAT(///5X, ’ [CS] Effective Stiffness Matrix’,//) 

8267 FORMAT (5X ,E10 . 3 , 3X ,E10 . 3 , 3X ,E10 . 3) 

8268 FORMAT (44X.E10. 3/ 

& 57X.E10.3/ 

& 70X.E10.3) 

END 


C 

C ############################################################################ 

SUBROUTINE BODNER (DSA , SA , STRC , NM , NS , INC , N , NN) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PURPOSE: Bodner-Partom Viscoplastic Model 

Note: 1) In This Subroutine, [SA] And [DSA] Contain The 

"Micro" Quantities For Aboudi’s Micromechanics Model 


2) Arrangement Of 

[DSA] & [SA] 

Variable 

Location 

1 

I Strain Rate 
1 

(1-6) 

1 

1 Stress Rate 

I 

(7-12) 

1 

1 Inelastic 
I Strain Rate 

i 

(13-18) 

1 

I 12 "Slots' 1 
I For State Variables 

i 

(19-30) 

1 

I Thermal Strain Rate 
1 ‘ 

(31-36) 

NM - MATERIAL NUMBER 
NS - SUBCELL NUMBER 


CALLED FROM: HYPELA (GMC) 
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c 

IMPLICIT REAL*8 (A-H.O-Z) 

C 

LOGICAL FORDIS .HEAT ,THME , NINDPT , PATRAN 
C 

COMMON /MLGIC/ FORDIS, HEAT, THME, NINDPT, PATRAN 

COMMON /MICMAT/ EA(2) ,ET(2) ,GA(2) ,FNA(2) ,FNT(2) ,ALPA(2) ,ALPT(2) , 
& VPBP(6,2) ,VPR0B(2,12) 

COMMON /MICRO/ VF , MATNUM (2 , 2) , NST , NMT , NB , NG , IDP , NSFD , NCMD 
COMMON /MICR02/ VF1 ,VC1 ,VF2,VC2 
COMMON /MICR03/ ICOUNT.NPLVL 
COMMON /MICR05/ CTEMP 
C 

DIMENSION SSC6) ,S(6) ,R(6) ,DSA(36) ,SA(36) ,STRC(6) 

C 



C Extract appropriate viscoplastic material constants 

C 

C 

DO = VPROB(NM.l) 

ZO = VPR0B(NM,2) 

Z1 = VPR0B(NM,3) 

BM = VPROB(NM,4) 

AN = VPROB(NM,5) 

Q = VPR0B(NM,6) 

C 



C Compute the deviatoric stress {S> in the subcell 



C 

TEMP = (SS(1) + SS(2) + SS(3))/3. 

S(l) = SS(1) - TEMP 

S(2) = SS(2) - TEMP 

S(3) = SS(3) - TEMP 

S(4) = SS(4) 

S(5) = SS(5) 

S(6) = SS(6) 

C 
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C - — 

C Predict the average plastic strain-rate 

C 

C 

AJ2=0.5*(S(1)**2+S(2)**2+S(3)**2)+S(4)**2+S(5)**2+S(6)**2 
SQ3AJ = SQRTC SS(1)**2 + SS(2)**2 + SS(3)**2 + 

& 2*(SS(4)**2+SS(5)**2+SS(6)**2) ) 

SQ2=1. 414215 
IF(SQ3AJ.EQ.O. ) THEN 
CALL ZER0R(R,6) 

ELSE 

R(l) = SS(1)/SQ3AJ 
R(2) = SS(2)/SQ3AJ 
R(3) = SS(3)/SQ3AJ 
R(4) = SQ2*SS(4)/SQ3AJ 
R(5) = SQ2*SS(5)/SQ3AJ 
R(6) = SQ2*SS(6)/SQ3AJ 
ENDIF 
C 

C - 

C If D0=0 then assume elastic and zero-out 
C [DSA(13-30)] (inelastic strain rate and 
C internal variable rates) , then return 

C 

C 

IF(DO.EQ.O) THEN 
DO 100 JJ=13,30 
DSA(JJ) = 0.0 
100 CONTINUE 

RETURN 

C - 

C Inelastic 

C 

ELSE 

C 

ZEF = ZO + q*SA(20) + 

& (1-Q)*(R(1)*SA(21)+R(2)*SA(22)+R(3)*SA(23)+ 

& R(4)*SA(24)+R(5)*SA(25)+R(6)*SA(26) ) 

C 
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IF(AJ2 .EQ. 0.) THEN 
AL=0 . 0 
ELSE 

ARG1=ZEF**2/ (3 . *A J2) 

C0N= .5*(AN+1 . )/AN 

ARG=C0N*(ARG1)**AN 

AL=D0/ (EXP (ARG) *SQRT(A J2) ) 

ENDIF 

C 

C 

C Inelastic strain rates 

C 

C 

DSA(13) = AL*S(1) 

DSA(14) = AL*S(2) 

DSA(IS) = AL*S(3) 

DSA(16) = 2*AL*S(4) 

DSA(17) = 2*AL*S(5) 

DSA(18) = 2*AL*S(6) 

C 

C 

C Plastic work rate 

c 

C 

WPD = S(1)*DSA(13) + S(2)*DSA(14) + S(3)*DSA(15) + 
& S(4)*DSA(16) + S(5)*DSA(17) + S(6)*DSA(18) 

C 

C 

C State variable rates 

C 

C 

DSAC19) = WPD 
Z0M=BM/Z0 

ZD=Z0M* (Zl-ZEF) *WPD 
DSA(20)=ZD 
C 

DSA(21)=ZD*R(1) 

DSA(22)=ZD*R(2) 

DSA (23) =ZD*R(3) 
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DSA(24)=ZD*R(4) 

DSA(25)=ZD*R(5) 

DSA(26)=ZD*R(6) 

END IF 
C 

RETURN 

END 

C ########################################################################## 

SUBROUTINE PLACE (A,B ,NN) 

C 

IMPLICIT REAL *8 (A-H.O-Z) 

DIMENSION A(NN) ,B(NN) 

C 

DO 10 1=1, NN 
B(I)=A(I) 

10 CONTINUE 
C 

RETURN 

END 

C ###################################################################### 

SUBROUTINE ZEROR(A.IDIM) 

C 

C PURPOSE: ZERO AN ARRAY 

C 

IMPLICIT REAL*8 (A-H.O-Z) 

C 

DIMENSION A(IDIM) 

DO 100 1=1 ,IDIM 
A(I) = 0.0 
100 CONTINUE 
RETURN 
END 

C################################################################## 
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